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1. Introduction 

Recent advancements of experimental techniques make it possible to fabricate 
nano-electronic devices where a quantum dot is connected to two superconducting 
electrodes. [T] Below the critical temperature, the electrons form a superconducting 
condensate (in other words, a single macroscopic quantum state). Therefore, in 
the case where the electrodes are superconducting, the quantum dot setup allows 
us to study the single electron tunneling between two condensates held at different 
chemical potentials, temperatures, or being forced to have different order parameters 
(e.g. different phases of the anomalous electron density). The mixture of different 
physical phenomena, such as single electron tunneling, quantum phase transition, and 
macroscopic condensation, opens the possibilities to study the fundamental physics. [2] 

The electron transport through a quantum dot involves three different energy 
scales: the tunneling coupling between the dot and the electrodes, the strength of 
electronic correlations inside the dot, and the order parameter for the superconducting 
state in the electrodes. Most of the theoretical research that has been done so far 
has been employing Keldysh nonequilibrium Green's functions (NEGF) or scattering 
theory type approaches. j3[ 21 [5] NEGF and scattering theory are able to treat the 
tunneling coupling exactly, but they usually fully neglect correlations inside the dot 
or they rely on mean field or perturbation theory to treat them. Here we develop 
an approach which is based on the Markovian quantum master equation [SJ [7] . The 
master equation approach to quantum transport works in the opposite regime - it can 
treat the correlations inside the dot very accurately (even exactly in the case of model 
systems) but the tunneling is usually considered in the Born-Markov approximation. 
Such an approach has been proved very useful for treating non-equlibrium transport 
problem in various quantum systems 0[l[ini[IIl[]]3[ni[Il[I]3[M 
been also applied to superconducting systems |20j . where the proximity effect in one 
dimensional wires was studied. A Lindblad master equation with quadratic Lindblad 
operators was obtained in the mean field approximation by mapping the many body 
super-operator to a single particle form [ST]. We note that a consistent treatment of 
the baths in the master equation approach poses many delicate issues [35] (see e.g. 
also discussion in Ref. |15j). 

Here we present a derivation of the master equation in the case when the 
electrodes are described by Bogoliubov-de Gennes Hamiltonians and then apply it to 
the non-equilibrium superconducting Anderson impurity model. We study in detail 
the transport properties of the model and the proximity effect in quantum dot. Three 
different regimes are considered. First, we focus on the generic case 2e + U ^ 0, 
where e is the resonance level energy and U the interaction strength, where an exact, 
analytic expression for the steady state is found. In the particle-hole symmetric 
regime 2e + U = we consider two cases, namely a dissipative one A < \e ± n/2\ 
and a non-dissipative one A > |e ± /j/2|, where /i is the chemical potential bias and 
A the magnitude of the superconducting order parameter. In the dissipative case 
the phase difference dependent non-equilibrium particle current, the energy current, 
and the proximity effect are obtained. In the non-dissipative case, the Josephson 
current originating from the Andreev bound states is discussed. The energies of the 
Andreev bound states and the corresponding particle current are obtained for arbitrary 
superconducting order parameter A and onsite energy level of the quantum dot e. 

The paper is organized as follows. In Sec. [2j we derive the master equation for 
a quantum system connected to superconducting baths, and then specialize on the 
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specific derivation for the out of equilibrium Anderson impurity model connected to 
the two superconducting leads. In Sec. |3j we present the numerical and analytical 
solutions of the master equation for the model cases. Conclusions are given in Sec. [4] 
We use natural dimensionless units throughout the paper, in which h = ks = |e| = 1, 
where — e is the electron charge. 



2. Markovian master equation for a quantum dot connected to 
superconducting baths 

In the derivation of the Lindblad master equation one usually assumes that the 
interaction operators between the system and the bath are written in a Hermitian 
form. [22] This is always possible and it usually simplifies the formal derivation. 
Therefore, we begin with the outline of a general derivation of the Lindblad master 
equation and highlight the main differences from the usual textbook approach. [22] 
The complete Hamiltonian is divided into three parts 

H = H S +H B +H h (1) 

where H-q denotes the bath Hamiltonian, H$ is the system Hamiltonian, and Hi is 
the interaction between the system and the bath. The interaction can always be 
represented in the following separable form 

Hi = Y, A ^B a (2) 

a 

where the operators A a (acting on the system) and B a (acting on the bath) commute 
[A a , B a ] — 0. As noted before, we shall avoid the common assumption that A\ = A a 
and fit = B a , since in our case the special form of the superconducting bath 
correlation functions induces two physically distinct contributions to the dissipator 
that are clearly separated only if we use the above form of the interaction 
However, since Hi needs to be Hermitian, the set {A a B a } has to include pairs of 
mutually Hermitian conjugate operators, i.e. for each a there exists a' , such that 
A a t — A^ al B a i = B^ a . The density matrix of the complete system satisfies the 
von Neumann equation. We use the standard Born-Markov approximations, namely 
that the denisty matrix of the complete system can be written in a separable form 
p{t) = ps <8> pe , where p§ denotes the density matrix of the system and p-Q the density 
matrix of the bath, which is assumed to be in a Gibbs state. Therefore, we can simplify 
the von Neumann equation and trace out the bath degrees of freedom. Further, 
by performing an additional secular approximation we obtain the Lindblad master 
equation for the reduced density matrix of the system 

-i[H hS ,p s (r)}+Vp s (r), (3) 

H hS =^2^2s a p(u)il- u (A a )fL(Ap) i (4) 

^Ps(r) = J2J2^^ (^UAp)ps(T)tl^(A a ) - {n_ w (^ Q )n w (^),p s (r)}), (5) 

where [•,•] denotes the commutator and {•,•} the anticommutator. The super- 
operators Ii u are projection super-operators on the eigenoperators of the system 
Hamiltonian H$ and are defined as 

fUO s )= Yl \e)(e\O s \e')(e% (6) 
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where |e)(e| are the projection operators on the possibly degenerate subspace of the 
system with the energy e (H s \e) = e|e)) and Os lS an arbitrary operator acting on the 
system. Note, that the equation ([3| can be brought to a standard Lindblad form since 
the matrix is Hermitian and positive semi-definite. 22J The functions S a p(uj) and 
7 a) g(w) are computed from the bath correlation function 

poo 

T a0 {uj)=r{B a ,B p \oj)= / dse™ s tr B (B a (T)B (T-s))=j a0 (uj) + iS a f ) (uj), (7) 

Jo 

la p{u) = 1 (B a ,Bp\u) = \{Y{B a ,Bp\u)+T*{Bl,Bt\w)), (8) 

S aP (u) = S(B a ,B p \oj) = i (r{B a ,B p \w) -r*(Bj,Bt| w )) , (9) 
where trs(») denotes a trace over the bath. 




Figure 1. Schematic illustration of the out equilibrium superconducting 
Anderson impurity model: one spin-degenerate level with energy e and local 
electronic repulsion U is connected to two supeconducting semi-infinite leads 
described by the Bogoliubov-de Genncs Hamiltonians. The chemical potential 
difference between left and right reservoir /i is included in the model through 
the bath correlation functions and describes the effect of a bias voltage, which 
is usually measured in experiment. In a similar manner the energy e may be 
associated with the gate voltage. 

Let us now apply the above consideration to a specific model, shown in Fig. 
[I] We consider a quantum dot connected to two uncorrelated one-dimensional 
superconconducting leads described by the Bogoliubov-de Gennes Hamiltonian 

h b=J2 e *( 6 £,At + 1>UaP-u) + H^b- k<i b ktt + e-'Xf&U,!)- (10) 

k 

Here &t _/i>fc,<r are creation/annihilation operators for an electron with spin a =t,i 
and single-particle energy e^, Ae 1 ^ is the complex order parameter, which governs the 
supeconducting properties of the leads. The index k runs over the modes of the left 
and right leads and A and phase <f> may have different values for the left and right leads 
(but we do not wish to burden the notation with additional indices). The quantum 
dot consists of one spin-degenerate level with on-site energy e and with local Coulomb 
interaction U > 0: 

H s = £^n a + Un-fUi, (11) 
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where 



= «t, 



is the number operator for electrons with spin a in the quantum 



dot. Here 4 and a a are creation and annihilation operators in the quantum dot, 
respectively. The interaction between the quantum dot and the superconducting leads 
is taken to be in the standard tunneling form 

Hi =J2 t '°A b l,<,a<T+ a ih,*)- (12) 



ka 



Since we are dealing with fermions the creation and the annihilation operators 
in the bath and in the system anticommute {a^ , frj^L} = 0. To establishe the 
connection with the master equation derived above ([3]), where we assumed that 
system and bath operators in Hi commute with each other, we perform a Jordan- 
Wigner rotation of fermonic creation and annihilation operators. Namely, we identify 
the operators in the interaction part of the Hamiltonian Hi as A a = a CT P B acting 
on the system and the corresponding B k , a = t k ^ a b\ CT P B acting on the bath, where 

Pb = exp (^m g b\ a b k ^ is the parity operator in the bath, which satisfies the 
following (anti)commutation relations 

[a(t),P B ]=0, {^,P B }=0. (13) 
It is easy to verify that A a and B k .a k (for k = 1, 2, . . . and er, a k =f , I) commute 



[A a ,B kt<Tk ] = [a a PB,t k ,a k bk,cr k PB} = tk,a k (ac,PBbk,a k PB ~ h,a k -PbOo-Pb) 



(14) 



= t k ,a k {acrPBbk,a k PB - a a PBb k ,a k PB) = 0. 

Hence, the interaction part of the Hamiltonian can be written as 

Hj = ((b{ a PB)(P B a a ) + (4P B )(P B ^ ff )) = J2(A*B k ^+AiBl a ). (15) 



fc.CT 



h . a 



Now we can calculate the correlation matrices Q for our model (see Appendix 



Appendix A). In order to obtain the dissipative part of the dynamics we have to 



find the projectors of the operators A a , A\ on the eigenoperators of the Hamiltonian 
as well. The eigenvectors (states) and the corresponding eigenvalues (energies) of the 
dot Hamiltonian are denoted as follows 



States 


|o> 


It) = 410) 


id = a }|o) 


|U) = a\a\ |0) 


Energies 





e 


e 


2e + U 



(16) 



Here the state |0) denotes the particle vacuum. Hence, the nonzero projections of the 
operators A a and A\ on the eigenspace of the Hamiltonian are 

n £ (a CT P B ) =\0){a\, tl e+u (a a P B ) = s a \a)(n\, (17) 

n_e(p B 4) = \<r){o\, n-e-t/(P B 4) = s«\n)(v\, 

were Sf = 1, sj, = — 1 and a denotes the opposite spin of a. Inserting the above 
projections (17) and the correlation functions calculated in the Appendix Appendix 
[A] into the master equation ([3]) we obtain the dissipative part of the Liouvillean of the 
quantum dot connected to a superconducting reservoir 



V (1 \p) = ]T ( 7 (1) (-c) {2\a){Q\p\G){a\ {|0)(0|,p}) 



(18) 



7 (l) (-e-£0 {mww){n-{\*m,p}) 

7 (1) (e) (2\0)(v\pW)(0\-{W)(vIp}) 
- 7 (1) (e + U) {2\a)(n\p\U){o\ {|tl)(tl|,p})' 
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and the Lamb shift term (see Ref. [22] for the definition) in the Hamiltonian 

#ls=E( SW(-e)\0)(Q\+SW(-e-U)\v)(a\ (19) 

+ S^(e)\a)(a\+S^(e + U)\n)(n\). 

In the particle-hole symmetric case (2e + U = 0) we have an additional contribution 
to the Lamb shift and the dissipator. This is a consequence of two effects: (i) 
the non- vanishing superconducting correlation functions r(PB&j. -p PB^_ k jjw) and 
r(6fc^Pe, &_fe 4.Pb|w), which signal a finite density of cooper-pairs and (ii) the twofold 
degeneracy of the energy zero in the dot, which ensures non-vanishing projections 
of the operators A a and A\ on the eigenspaces of the Hamiltonian H$ with 
opposite energies. Therefore, products n e (a cr PB)n_ c (a ff PB) and n e (PBc4)IL c (PBc4) 
appearing in the sums Q and ([5| do not vanish as in the non-degenerate case 
(2e + U 0) and we obtain the following, additional contributions to the dissipator 

V^(p) = J2( 2 7 (2) ( e )|U}Hpk><0| (20) 

+1 {2 \-e) {2\a)mnm-{\n){%p}) 

+ 2 7 ( 2 >*( e )|0)(a|p|a)(t|| 
+ j^*(-e) (2\*){U\p\0}(cr\-{\a}{n\,p}) 

and to the Lamb shift 



< = 2 



i(5( 2 )(-6)|n><0| + ^ 2 )*(-e)|0)(tl|). (21) 



For the sake of simplicity the above expressions for dissipators ( 18 ) , ( 20 1 and the Lamb 



shifts (21), (19) are written for one bath only. The contribution of the second bath is 
identical and additive, so the total dissipator and the Lamb shift become 

V = V L +V R , H LS = Pls.l + Pls,r, (22) 

where v = L, R. Note that we neglect the broadening of the systems energy levels 
due to the coupling to the leads, i.e. the levels are infinitely narrow. The broadening 



can be included by hand setting for example rj — k (see Appendix Appendix A ) or by 
self-consistent treatment of the master equation as suggested in Ref. [23] ! 

3. Solution of the master equation 

In this section we shall find the steady state density matrix pness of the master 
equation (J3j) . First we consider the non-degenerate quantum dot, where we have 



only one contribution to the dissipator, namely (18), and the Liouville equation is 
simplified to a rate equation. An explicit analytic form of steady state is obtained. In 
the second subsection we consider the particle-hole symmetric case, where the steady 
state is calculated numerically. We find nontrivial non-equilibrium sub-gap dynamics 



due to the effect of the Lamb shift ( 21 ). In both cases we discuss the non-equilibrium 
particle current and energy current defined as a change of the number of particles in 
the system and system's energy, respectively, due to the interaction with the left bath 

J n = t>l{n) +i[# LSl L, n], n = J2 a W, (23) 
J e =T$(H S ) +i[H hS<u H s ], (24) 
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where the superscript H denotes the Heisenberg representation of the superoperator 
T>^. We also discuss the proximity effect, namely the cooper pair density in the 
quantum dot 



A dot e^ 



(ar a ;>- 



(25) 



3.1. Non-dengerate quantum dot: 2e + U 7^ 



As already explained, in the non-degen erat e case we need to take into account only the 
first part of the dissipator (2?W), Eq. (18). The subindex v = L,R in the correlation 



(i) (i) 

functions 7,7 and Si, denotes different baths. In this case the steady state can be 
found analytically by writing the Liouvillean in the matrix form and noting that the 
coherences decouple from the rates, which results in a simple rate equation the solution 
of which is 

Pness 
Po 

Pi 

P2 



(po|0}<0| +pi(|t><t| + + P2 \U)(n\)/(po + 2pi +P2), (26) 

+ 7 R 1) «)) (ll 1} (U + e)+ >yg>(U + e)) , 

(-e) + 7^ ("*)) (7^ (U + e) + 7 £> (U + ej), 



Interesting observables in the non-degenerate case are the particle current (23) and 
the energy current (|24|), which simplify to 



J" 



E (-47[ 1) (-e)|0)(0| + (-2 7 [ 1) (- 
+ 4^\e + U)\n){U\] 



and 



E (-4 £7 [ 1) (-e)|0)(0| + (-2[/7[ 1) (^ 

a 

+4(U + e) 1 i 1 \e + U)\n)(n\] 



U) + 2^\e))\a)(a\ (27) 



U) + 2 n [ 1) (e))\a)(a\ (28) 



respectively, where the subscript (L) denotes the left bath. By using the equations 
(26 1, (p7| , and ( 28 ) we can easily calculate the expectation values of the particle and 



energy currents in the steady state. The chemical potential is included by replacing 
£^e±/i/2in the left (+) and the right (-) bath correlation functions. The qualitative 
behavior of the currents (and the differential conductance) can entirely be explained by 
the electronic density of states in the superconducting leads (superconducting density 
of states - SDOS), 

PL H = e(M-A)M/^ 2 -A2, (29) 

where O(w) is the Heaviside step function (shown in Fig. [I]). The couplin g strength 



to the baths 7^^ is proportional to the SDOS, as shown in the Appendix Appendix 
|A| . Hence, the main features of SDOS, namely a gap 2A where the electronic 
density of states is zero and a divergence at the border of the gap are reflected 
in the current-voltage characteristics (see Fig. [2| and the differential conductance 
G = d(J n )/dii map (see Fig. [3]) . Far from the gap, i.e. when Al,r -C |e ± p/2\ and 
Al,r <C |e ± fJ,/2 + U\, the current approaches the value calculated for the normal 
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leads. As we approach the superconducting gap (by changing the chemical potential 
fi or the onsite energy e) we observe a peak in the differential conductance when one 
set of the following conditions is satisfied 

u = A L - /x/2 and uj < -A R + fi/2 or u> > A L - fi/2 and uj = -A R + fi/2, (30) 

where w = eorw = e + J7is the transition frequency between the subsequent levels of 



the dot. The conditions (30) are valid if /i > 0, but for fi < the roles of the baths are 
exchanged (Al <-> Ar). After the peak we observe negative differential conductance 
as a consequence of the decreasing density of states, which is clearly shown in Fig. 
[2] The characteristic distances appearing in the differential conductance map in Fig. 
p3] can easily be calculated from the first and the last condition in ( 30 ) , namely the 



interaction energy U, the sum of the superconducting order parameters in left and 
right lead Al + Ar , and the difference of the superconducting order parameters in 
the leads |Al — A R |. The interaction energy determines the difference between the 
two possible transition energies U) for one particle transfer and thus also the relative 
shift of the diamond structures on the e (gate voltage) axis in Fig. [3j The distance 
Al + Ar determines the bias voltage (or chemical potential) that has to be applied 
to the leads in order to get the maximal particle current, namely to align the top 
of one superconducting gap with the bottom of the other; they align at the energy 
|Al — Ar|/2, which is shown in panel b) of Fig. [2j 

Note that, the cooper pair density in the quantum dot (the proximity effect) 
vanishes since there are no coherences in the steady state of the non-degenerate dot. 
If the level is placed inside of the superconducting gap of both superconductors the 
dot is not coupled to the environment and we obtain a trivial unitary evolution. In the 
next subsection we shall show that in the particle hole symmetric case the evolution 
of the level inside of the gaps is changed due to an additional non-vanishing Lamb 
shift term. 



a) 
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M/2 
f/2 



\n) 

.JO) 

:(A L +A R )/2 



:]= |A L -A R |/2 



Figure 2. a) Current- voltage characteristics of the Anderson impurity. Dashed 
line corresponds to the normal leads (Al r = 0) and the full line corresponds 
to the superconducting leads (Al,r = 0.5). For superconducting leads the shift 
of the step by 2A and the negative differential conductance are a consequence 
of the SDOS. Other model parameters: e = 1, U = 2,T l ,r = 0.1, T = 0.1. b) 
Schematic representation of the configuration with the maximal current. The 
chemical potential is /i = Al + Ar, therefore the top of right superconducting 
gap aligns with the bottom of the left superconducting gap; they align at the 
energy e = |A L - A R |/2. 
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Figure 3. Density plot of differential conductance (G = d{J n ) /d/i) as a function 
of the level position e and the bias voltage (chemical potential difference) fi, for 
different values of the order parameters Al.r, (indicated on top of each panel). 
The characteristic distances denoted in the first three panels are calculated from 
the equation J30l; see main text and panel b) of Fig. [2] Dark/light gray represent 
high/low values of differential conductance G with, in each plot, suitably adjusted 
relative scale. Other model parameters: U = 2, Tl,r = 0.1, T = 0.1. 



3.2. Particle-hole symmetric case 2e + U = 

In the symmetric case (2e + U = 0) an additional nontrivial term appears in the 
Lindblad master equation, which is not present for non-superconducting leads and 
describes Cooper pair tunneling between the bath and the dot. Not surprising, the 
dependence of the time evolution on the phase of superconducting order parameters 
in the leads is introduced through this second part of the dissipator ( 20 ) . An analytic 



solution in this case is cumbersome, since the populations are now coupled to the 
coherences through the extra dissipator, therefore we rely on numerical calculations 
to find the exact steady state. As in the non-degenerate case we study the particle 
current, energy current and the proximity effect. It is interesting, that in equilibrium 
we have no particle current although the superconducting parameters in the left 
and the right lead have a different phase. The energy current and the Cooper pair 
density in the dot are zero as well. However, in out of equilibrium steady state 
the currents depend on the difference of the superconducting phases in the leads 
Acf) = 4>h — 4>R ( see Fig. [4]) . Moreover, a finite non-equilibrium proximity effect ( 25 1 
is obtained; also shown in Fig. [4] Although we are unable to derive exact analytic 
expressions for the observed quantities we find some interesting effects when passing 
form the small bias regime /i < 2(|e| + A), where Al = Ar = A, to the large 
bias voltage regime /i > 2(|e| + A). If wc change the chemical potential difference 
from small to large values we observe a phase flip, which is it— shift in the phase 
dependence of the order parameter in the quantum dot. For zero bias voltage the 
phase profile is linear and it is antisymmetric under the reflection around the point 
(Aci, </>dot) = ( 7T , 7r )- This symmetry in the phase dependence 0dot(A0) is present 
whenever the temperature difference or the chemical potential difference between the 
baths is zero. This is illustrated in Fig. [5j where we show the dependence of the phase 
0dot and the amplitude Adot of the complex order parameter in the quantum dot on 
the superconducting phase difference Acf>. 

Next we consider the non-dissipative case, namely a narrow dot level (i.e. when 
the life time of the electron on the level is much larger than 1/2A) is placed inside 
the superconducting gap below the chemical potential of the baths. Surprisingly, we 
observe a finite non-dissipative particle current, which comes from the distortion of the 
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Figure 4. Out of equilibrium phase difference A0 = 0l — 0r dependent 
particle current (J n ) (dashed line), energy current (J e ) (dash dotted line) and 
proximity effect (Ad ot = |(a^a^}|, full line). The left panel is calculated for a 
temperature bias (fi = 0, Tl = 0.2, Tr = 1) and the right panel shows the results 
obtained for a nonzero bias voltage (or chemical potential; /i = 4, Tl.r, = 0.2). 
Note, that if Tl = Tr and fi = the particle current, energy current, and 
the Cooper pair density in the quantum dot vanish, therefore the obtained 
phase dependence is a purely non-equilibrium effect. Other model parameters: 
2t = -U = -2, T = 0.1, A l ,r = 0.5. 




Figure 5. The proximity effect in quantum dot A ( j ot e"' >dot = (a-j-a^) for various 
regimes. The color and thickness distinguishes between the temperatures of the 
right bath (thin, black) Tr = 1, (thick, orange) Tr = 0.1 and the dashing 
represents the chemical potential (full line) /i = 0, (dashed lines) fi = 0.5, (dotted 
lines) /i = 4. The curves calculated for large chemical potential bias are shifted by 
7r relative to the curves calculated for small chemical potential bias. An additional 
small phase shift occurs if a chemical potential bias and a temperature bias are 
present. Other model parameters: e = — 1,T = 0.1, Tl = 1, Al,r = 0.5. 



Hamiltonian due to the coupling to the baths i.e. from the Lamb shift. The modified 
Hamiltonian in the described case is 



( 2SW(- 



V (- 





2S«(e) 








25( 1 )(e) 




\ 






2S«(- e ) j 



(31) 



SM(U) = S¥ ) (") + Sg\u), j = 1,2 



The Hamiltonian in the degenerate subspace spanned by the states |0) and | tl) is 
perturbed by the interaction with the environment. This lifts the degeneracy and 
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obtained energy eigenstates 
|* ± ) = -L(± e -^/ 2 m) + |0)), E ± =2S^(-e)±4\S^(-e)\ C os(^\ (32) 

are no longer eigenstates of the particle number operator. The particle current in this 
case simplifies to 



J" = 4i (-4 2) (-e)|n)<0| + sW*(-e)\0)(U\) ■ 
and the expectation value of the particle current in the states \^±) is 

(* ± |J n |*±)=T2|5( 2 )(- e )|sin^V 



(33) 



(34) 



In the above equations (32 34J) we assume that Al = Ar = A and |e| = U/2 < A. 
The current carying states can be interpreted as the Andreev bound states (see Ref. 
[5]). Further, we assume that the dot is in the Gibbs state pc of the modified 
system Hamiltonian H sym = H$ + H^s and calculate the particle current (see Fig. [6| . 
Interestingly, the particle current oscillates around zero also in the case where we have 
a non-zero chemical potential difference. 





1 2 3 4 5 6 

\<!> 

Figure 6. We plot current carrying energy levels of the Lamb shift modified 
Hamiltonian H sym (left) and the corresponding particle current (right) calculated 
in the Gibbs state of H Byul . Black dashed lines show the equilibrium result, namely 
the analytic result of equations i \'V2j and | |34| , calculated for Al,r = 2, fi = 0. Out 
of equilibrium (orange lines) crossing of Andreev bound state energies is avoided 
(parameters: Al = 2, Ar = 5,/i = 1). Note, that in case of nonzero chemical 
potential difference or bias voltage (orange lines) the particle current oscillates 
around zero. Other model parameters: e = — 1/2, T = 0.1, Tl r = 1. 



4. Conclusions 

We derived a master equation for the electron transport through the quantum dot 
connected to two superconducting leads. In our derivations the Born-Markov and 
the rotating-wave approximations were used, which reduce the master equation to the 
standard Lidblad form with nontrivial dissipators and the Lamb shift terms originating 
from the superconducting baths. Then, the master equation was explicitly solved for 
out of equilibrium Anderson impurity model and the exact steady state density matrix 
in the generic regime 2e + U =/= was found. In the particle-hole symmetric regime 
2e + U — a phase dependent dissipator was obtained. Surprisingly, the equilibrium 
solution in this case does not exhibit a superconducting phase difference dependent 
particle current, whereas in the out of equilibrium steady state the particle current, 
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the energy current, and the Cooper pair density of states in the quantum dot depend 
on the difference of the superconducting phases in the baths A</>. In the sub-gap 
case, Andreev bound states [22J are found as eigenstates of the Lamb shift perturbed 
Hamiltonian. Their energies and the corresponding non-dissipative particle current 
were obtained also in the non-equilibrium situation. The master equation derived in 
this article can be extended to treat larger systems, e.g. the double quantum dot, 
molecules, and one dimensional wires. 

This work has been supported by the Francqui Foundation, Programme d'Actions 
de Recherche Concertee de la Communaute francaise (Belgium) under project 
"Theoretical and experimental approaches to surface reactions", the grants P 1-0044 
and Jl-2208 of Slovenian Research Agency (ARRS), and a scholarship from the Slovene 
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Appendix A. Calculation of correlation function for supeconducting bath 

In the appendix we shall obtain the correlation function in the bath with the 
Bogolyubov-de Gennes Hamiltonian, which can be diagonalized by the Bogolyubov 
transformation 

bfc.t = - u(e k )d kA + v(e k )d\ x , 6_ fc4 = u(e k )d k ,i + v(e k )d\^ (A.l) 



2 2u(e)' w y 2 2cu(e) 
to(e k ) = Je 2 k +A 2 , tan(20 fc ) = - — . 
H B = ^y^u(e k )dl a d K<T . 

k a 

We assume that the bath is in the equilibrium at temperature T 

(4,<Av>b = «(efc)4,fc"W, n(e k ) = l + ^ l{tk)/T - (A.2) 
Hence, the only nonzero correlation functions are 

IWLA^bM =iK £fc )| 2 "fr^ ■ +ib(^)| 2 1 ~ n( ; fc ], , (A.3) 

uj + uj(e k ) + in ui - LO{e k ) + 177 

%/B,p B ik) ^ii^ £fc )i 2 "? £fc L. +^fa)i 2 1 ' n( ; fc ], , 

w + cj(e fc ) +177 w - co(e k ) + irj 
T(b k , t P B ,b. kii P B \u) = -iu(e k )v(e k ) ( " (cfc) , ) , 

r(p B b{ v P B bi k ^) = iu(e k yv(e k y ( 1 7 ( ; fc ? ■ ) • 

V^ + ^l 6 *;) + 17 7 oj - ui{e k ) +in J 
At the end of the calculation we shall take the limit i] — > + . Further, we assume an 
energy and spin independent coupling to environments k — nJ2k \tk,a\ 2 S(e — e k ) and 
obtain 

7 (1) (^) = ]>> fe , ff | 2 7 (ft< CT A,.PBM 

= £ [°° de6(e-e k )\t k ,^( ll«(^)l 2 »* , ^)| 2 (l-« fe 
,. J— 00 



77 2 + (wfc + w) 2 7] 2 + (w - w fe ) 



W-oo \V 



\u(<)\ 2 ni<) | Ke)| 2 (l-n(e)) - . y 4) 



■ 2 + (w(e) + w) 2 r/ 2 + (u;-u;(e)) 2 

K f°° ( nit) 1 - n(e) 

1 den ' 



7T 



ry 2 + (uj(e) + uj) 2 i] 2 + (oj — to(e)) 



7 (1) M = lim 7 (1) (^, V) = Kp L (u>)(l - n(u>)), 

rj— >0 

7 (2) (w,r?) = ^|^| 2 7(6 fc , t P B ,fo-fe 4 PBk) 
fc 

k f°° ( n(e) 1 - n(e) 

«-/ deu e)7j(e )n - --- 

Wo \i] 2 + (oj(e) +uj) 2 n 2 + {uj - u(e)) ■ 

7 (2) (w) = lim 7 (2) (w,r/) = 2k Pl {uj){1 - n(u)))u{u)v{u))(Q{u)) - 6(-w)), 
77— >o 
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where pl{u) = 0(w — A)w/\/w 2 — A 2 is the superconducting density of states and 
8(w) is the Heaviside step function. All other correlation functions are up to a sign 
equal to ^^(uj) or ^^(lu). In order to determine the Lamb shift we have to calculate 
the following sums 

S^\oj) = J2\ t kM 2 S(PBbl a ,b Ka P B \u) = Y,\ t kM 2 S(bk,aPB,PBb{joj), (A.5) 
fc k 

S&(u;)=J2\ t k,a\' 2 S(b k , t P B) b- k , l P B \u J ), 
k 

k 

J2MS(P B bl >v P B bl kti \u) = -S^(lu), 

k 

E Ka\ 2 S(P B bl v P B bl ktt \ij) = SW*(U>). 
k 

This can be done numerically using the relations (A.3) and the definitions Q. The 
results are independent of the used bandwidth in the sums (A.5). 



